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Abstract 

In usual stochastic volatility models, the process driving the volatility of the asset price evolves accord- 
ing to an autonomous one-dimensional stochastic differential equation. We assume that the coefficients 
of this equation are smooth. Using Ito's formula, we get rid, in the asset price dynamics, of the stochastic 
integral with respect to the Brownian motion driving this SDE. Taking advantage of this structure, we 
' propose 

- a scheme, based on the Milstein discretization of this SDE, with order one of weak trajectorial 
convergence for the asset price, 

- a scheme, based on the Ninomiya-Victoir discretization of this SDE, with order two of weak con- 
vergence for the asset price. 

We also propose a specific scheme with improved convergence properties when the volatility of the asset 
price is driven by an Ornstein-Uhlenbeck process. We confirm the theoretical rates of convergence by 
Qh ' numerical experiments and show that our schemes are well adapted to the multilevel Monte Carlo method 

^ (-h introduced by Giles [2008 b, 2008 a]. 
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Introduction 

On ■ 

There exists an extensive literature on numerical integration schemes for stochastic differential equations. 
To start with, we mention, among many others, the work of Talay and Tubaro [1990] who first established 
an expansion of the weak error of the Euler scheme for polynomially growing functions allowing for the use 
of Romberg extrapolation. Bally and Taley [1996] extended this result to bounded measurable functions and 
Guyon [2006] extended it to tempered stable distributions. More recently, many discretization schemes of 
higher weak convergence order have appeared in the literature. Among others, we cite the work of Kusuoka 
[2001, 2004], the Ninomiya and Victoir [2008] scheme which we will use hereafter, the Ninomiya and Ninomiya 
[2009] scheme and the scheme based on cubature on Wiener spaces of Lyons and Victoir [2004] . 
Concerning strong approximation, the Milstein scheme has order one of strong convergence. Unfortunately, 
it involves the simulation of iterated Brownian integrals unless a restrictive commutativity condition is 
satisfied. Under ellipticity, Cruzeiro et al. [2004] have recently proposed a discretization scheme which gets 
rid of these iterated integrals and has nice strong convergence properties. More precisely, for each number of 
time steps, there exists a Brownian motion different from the one giving the Brownian increments involved 
in the scheme such that the strong error between the scheme and the stochastic differential equation driven 
by this new Brownian motion is of order one. We call such a property weak trajectorial convergence of order 
one. Weak trajectorial error estimation is exactly what is needed to control the discretization bias for the 
computation of path dependent option prices. 
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Stochastic volatility models, which have now become a standard of the market, are an eloquent example of 
the use of stochastic differential equations in finance. In our study, we will consider the following specification 
of a stochastic volatility model for an asset (S^tgfo.T] : 



dS t = rS t dt + f(Y t )St(pdWt + y/T=?dB t y, S*o-s o >0 (1) 
dY t = b(Y t )dt + a(Y t )dW t ; Y = y a 

where r the instantaneous interest rate, (-Bt)te[o,ri an d (W / t)te[o,T] are independent standard one-dimensional 
Brownian motions, p G [—1,1] is the correlation between the Brownian motions respectively driving the asset 
price and the process (it)tg[o,Tl which solves a one-dimensional autonomous stochastic differential equation. 
The volatility process is (/(5 / t))tG[o,T] where the transformation function / is usually taken positive and 
strictly monotonic in order to ensure that the effective correlation between the stock price and the volatility 
keeps the same sign (the function a usually takes nonnegative values). This specification nests almost all 
the known stochastic volatility models: 

• Hull&Whitc [1987] model (p = 0) and Wiggins [1987] (p ^ 0) 

dS t = rStdt + y/ZSt^pdWt + y/T^dB^ 
dY t = pY t dt + (Y t dW t 

which can be expressed as (JTJ with f(y) — ^/y, b(y) = py and a(y) — (y. Nlote that it can also be 
seen as (JTJ> with f(y) = e y , b(y) = § - ^ and a{y) = §. 

• Scott's [1987] model which generalizes Hull&White model 

dS t = rS t dt + a e Y <S t ( K pdW t + ^T~fdB^ 

dY t = k(9 -Y t )dt + vdW t (2) 
f(v) = ', b(y) = k(6 - y) and a(y) = v. 

• Stein&Stein model [1991] 

dS t = 
dY t = 

• Quadratic Gaussian model 

f dS t = 

1 dY t = 

• Heston's [1993] model 

f dS t = rStdt + JTtSt^pdWt + yJT^fdBt) 

\ dY t = K{6-Y t )dt + v^/YtdWt 

f(y) = Vv> b (y) = K (Q _ v) and a (y) = u Vy- 



rStdt + Y t S t [pdWt + ^T~i?dBt) 
k(6 - Y t )dt + vdW t 

= y, b (y) = K (o - y) and <r(y) = v - 



rStdt + Y t 2 S t (pdWt + ^/T~7dB t j 
k{6 - Y t )dt + vdW t 

= y 2 , b (y) = K (Q - y) and <r(y) = v - 
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In all but the last example, the volatility of the asset is driven by an Ornstein Uhlenbeck process. 

The development of specific discretization schemes for stochastic volatility models has only received little 
attention. We mention nevertheless the work of Kahl and Jackel [2006] who discussed various numerical 
integration methods and proposed a simple scheme with order 1/2 of strong convergence like the standard 
Euler scheme but with a smaller multiplicative constant. Also the numerical integration of the CIR process 
and of the Heston model received a particular attention because of the inadequacy of the Euler scheme due 
to the fact that both / and a are equal to the square root function (see for example Deelstra and Delbaen 
[1998], Alfonsi [2005], Kahl and Schurz [2006], Andersen [2007], Berkaoui et at [2008], Ninomiya and Victoir 
[2008], Lord et al. [2008], Alfonsi [2009]). An exact simulation technique for the Heston model was also 
proposed by Broadie and Kaya [2006]. 

In the present paper, we assume in return that the functions /, a and b are smooth and do not deal with 
the Heston model. Our aim is to take advantage of the structure of ([TJ to construct and analyse simple and 
robust ad'hoc discretization schemes which have nice convergence properties. 

For a start, we make a logarithmic change of variables for the asset : the two-dimensional process 
(X t :— log (S t ) , Yt)te[o,T] solves the following SDE 

dX t = ( r -±f(Y t ))dt + f(Y t )(pdW t + ^/r^dB t y, X = \og( So ). (g) 
dY t = b(Y t )dt + a{Y t )dWu Y = y 

Our main idea is to get rid in the first equality of the stochastic integral involving the common Brownian 
motion (Wt)te[o,T]- In a U what follows, we assume that 

(A) f and a are C 1 functions and a > 0. 

One can then define the primitive F(y) = £(z)dz and apply Ito's formula to get 

dF(Y t ) = -(Yt)dYt + \{af - fa')(Y t )dt. 
a Z 

Therefore (X t , F t ) te [o,T] solves 

dX t = pdF(Y t ) + h(Y t )dt+^T~p^f(Yt)dB t 
dY t = b(Y t )dt + a{Y t )dW t 

where h : y H> r — \f 2 {y) — p(-f + h( a f ~ f a '))(y)- We discretize the autonomous SDE satisfied by Y 
using a scheme with high order of strong or weak convergence depending on wether one is interested in 
path-dependent or vanilla options. Then, in the dynamics of X, we only need to discretize the standard 
integral h(Y s )ds and the stochastic integral f(Y t )dB t where (Yt)t£[o,T] an d (£?t)tg[o,T]- 

We recall that usual weak convergence is the right notion to analyse the discretization bias for plain 
vanilla options whereas weak trajectorial convergence permits to deal with path-dependent options. The 
first section of the paper is devoted to path-dependent options. Combining the Milstein discretization of the 
one-dimensional SDE satisfied by (Xt)te[o,T] with an appropriate discretization of the integral J Q T f(Y t )dB t 
based on the independence of (^t)te[o.T] an d (£?t)te[o,T]: we obtain a scheme with order one of weak tra- 
jectorial convergence. In the second section, using the Ninomiya- Victoir discretization of the SDE satisfied 
by 0"t)te[o,Tl) we construct a scheme with order two of weak convergence. Since the SDE satisfied by Y is 
one-dimensional, the Ninomiya- Victoir scheme only involves two one-dimensional ODEs whose solutions are 
available in closed form. The last section is devoted to numerical experiments which confirm the theoretical 
rates of convergence. We also show that our schemes are well adapted to the multilevel Monte Carlo method 
introduced by Giles [20086, 2008 a]. 
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Notations 

We will consider, for a number of time steps N > 1, the uniform subdivision J\ N = {0 = to < ti < ■ ■ ■ < 
tjy = T} of [0,T] with the discretization step 5n = -j^- 

We denote by ip the greatest lower bound of the function ifi : y H> f 2 (y) and by tp its lowest upper bound. 
We also introduce the following notation : 

-0(y) = 




1 An efficient scheme for path dependent options pricing 

Building a first order strong convergence scheme for a two dimensional SDE is not an obvious task. Even 
the ad'hoc schemes provided by Kahl and Jackel [2006] exhibit a strong convergence of order ^ . 

Actually, the natural candidate for this purpose is the Milstein scheme. Unfortunately, the commutativity 
condition which permits to implement it amounts to of = in our setting. This condition is typically true 
when either / is constant or o = 0. Both cases are of no practical interest since they lead to a deterministic 
volatility. 

However, since the inherent Brownian motion is not essential for applications in finance, the usual strong 
convergence criterion is not adapted for estimating the error of a scheme in pricing a path dependent option. 
What is more relevant is the approximation in law of the whole trajectory of the process considered for 
instance by Cruzeiro et al. [2004] . Using an ingenious rotation of the Brownian motion, these authors have 
constructed a discretization scheme allowing for a weak convergence on the whole trajectory of order one 
which avoids the simulation of the iterated stochastic integrals. 

For the SDE ©, the discretization scheme of Cruzeiro, Malliavin and Thalmaier writes as 

Xf™ = Xg MT +(r- H^fll^ 5n + pf(Y™ T )AW k+1 + %of (Y™ T )AWt +1 

+^7of'(Y t c h MT )AW k+1 AB k+1 + VT^7f(Y t c k MT )AB k+1 - ^of(Y t c k MT )AB 2 k+1 

(5) 

= Y t c k MT + (b{Yf k MT ) + \(^f - oo')(Y t ^)) S N + o(Y™*)AW k+1 

+ \oo'{Y t ^)AWl +l - ^fABl +1 

where AW tk+1 — W tk+1 — W tk and AB k+ i = B tk+1 — B tk correspond to the Brownian increments. 

We set out to construct a much simpler scheme having the same order of weak trajectorial convergence 
by taking advantage of the particular structure of the SDE defining stochastic volatility models. We first 
begin with the general case of any process (^t)te[o.T] driving the volatility and then consider the case of an 
Ornstein-Uhlenbeck process where we obtain more precise results. 



1.1 General case 

A discretization scheme will naturally involve the Brownian increments. Thanks to the independence 
between (lt)te[o,T] and (-Bt)te[o,T]j we can construct a vector (X to , ■ ■ ■ , X tN ) using only (ABi , . . . , ABn) 
and (Yt) t £[o,T]j which has exactly the same law as (X to , . . . ,X tN ) : 



Lemma 1 — V0 < I < N, let v t = ^ f^ +1 ip(Y s )ds. The vector (X t „, . . . , X tN ) defined by 

X t„ = X tg 

pt k fc-1 

Vl<k<N 1 X tk =X t0 +p{F{Y tk )-F{Y t0 ))+ / h{Y s )ds + ^/l - V^T AB Z+ 
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has the same law as (X to , . . . , X tN ) . 

Proof : The proof is elementary. Conditionally on Y, the two vectors are Gaussian vectors with the same 
mean and covariance matrix. □ 



In order to approximate (X tk )o<k<N , one needs to discretize Vk for fc G {0, . . . , N — 1}. If (v k )o<k<N-i 
is an approximation of (vk)o<k<N-i, then by Doob's inequality 



sup 

0<fc<JV-l 



E 

\Z=0 



N-l 



fc=0 

1 Ar_1 

^ ^E E [(^- y . 



fc=0 



as soon as ip = ini x ip(x) is assumed to be positive and, V0 < k < N — 1, is greater than Consequently, 
to obtain a scheme with order one of strong convergence for (X tk )o<k<N , one needs that V0 < k < N — 
= O ( a75-)- According to the treatment of the term I J 2 defined by ([9|) in the proof of the 



1,E 



Theorem [1 below, one has V0 < k < N - 1, 



v k -(iP(Y tk ) + 



>N 



(W s -W tk )ds 



= o 



N 2 



(6) 



This equality still holds true when replacing Y by a scheme with order one of strong convergence in the 
term with sign minus of the left hand side. Better still, (^F(Y tk ) + h(Y s )ds^j is approximated with 

strong order one when replacing Y by such a scheme and using a rectangular discretization for the integral 
in time. 

For all these reasons, we choose the Milstein scheme for Y : 



VQ < k < N - 1, Y t N k+i = Y» + b(Y£)6 N + o-(Y»)AW k+1 + -aa'{Y t N h ) (AW 2 k+1 S N ) 
and we write our scheme as follows 



WeakTraj_l scheme 



K +1 = K + P [ nYf k+1 ) - F(Y« ) ) + 6 N hQ% ) 



N- 



vr 



A 



~ Ar a^'(Y N ) /■**+! \ 
1>(X£) + 6n " J (W s - W tk )dsU± AB k+1 



t =vo- 



(7) 



Note that in order to implement this scheme, one needs to simulate both the Brownian increment AWk+i 
and the random variable J. fc+1 (W s — Wt k )ds. This is straightforward as one can easily check that 



V Jt 



AW k+1 
(W s -W th )ds 



)■"(( 



Sn S%/2 
S 2 N /2 4/3 



We can now state our first main result : 
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Theorem 2 — Under the assumptions of Lemma\3\ and if 
(HI) f and a are C 3 functions, — and f f are bounded 
(H2) V>>0 

(H3) there exists a constant K\ such that, V(x, y) € R 2 , 



{bh' + -h"){y) <K x {l + \y\) 
<rh'(y)\ <Xi(l + |y|) 
h(y)-h(x) <Ki\y-x\ 



(H4) there exists a constant K2 such that, y(x,y) € R 2 , 

(b^' + ^")(y)\<K 2 (i + \ y \) 

a^'(y) - o-ip'(x)\ < K 2 \y - x\ 

then the WeakTraj-1 scheme has order one of weak trajectorial convergence. More precisely, for each p > 1, 
there exists a constant C independent of the number of time steps N such that 



E 



max 

0<k<N 



(x tk ,Y th ) - (X»,Y» 



2p 



< 



c 



The proof of the theorem relies on the order one of strong convergence of the Milstein scheme (see Milstein 
[1995] for the particular case p = 1) : 

Lemma 3 — Suppose that 

(H5) b and a are C 2 functions with bounded first and second derivatives 
(H6) there exists a positive constant K such that V(x, y) € R 2 

\ao / (x)-aa / (y)\ < K\x - y\ 

then, Vp > 1, there exists a positive constant C p independent of N such that 

( max Y u . - % N 2P ) < C p 6%. 



E max 

\0<k<N 



The proof for general p is postponed to the appendix. 

Remark 4 — Before giving the proof of the theorem, we make a few comments on its assumptions. (Hlty 
implies that h and ip are C 2 functions which was implicitly assumed in (T^) and (H\Q). The latter assumptions 
are expressed in a reduced form. One can check that the following conditions on the coefficients of the original 
SDE are sufficient for them to hold : 

• / and a are bounded C A functions with bounded derivatives. 



b is a bounded C 3 function with bounded derivatives. 



G 



• 3<7o > such that Vy € R, a(y) > ctq. 

Proof of the theorem : Throughout the proof, we denote by C a constant which can change from one 
line to another while always being independent of N. Thanks to Lemma[3j we just have to control the error 
on X : 



m&x\X tk -X t N \ 2p 

0<k<N k 



where 



and 



max 

0<k<N 


k-l 1 

p(F(Y tk )-F(Y»))+J2( 


f J+1 h{Y s )ds - S N h(Y f f) 




1 ' 3+1 tP(Y s )dsAB j+1 




5n J 





\ 



' ~ aip' (Y t N ) /■*<+! \ 



m j tj 



< 3 2 P -l (p 2 P/o+/i + (1 _ p 2 )P/2) 



io = E 



max 

0<fc<JV 



F{Y tk )-F(Y t 



2p 



A = E 



I E ( f 1+1 h(Y s )d S - S N h(Y t f) 

3=0 \ Jt i 



2p 



E 



max . , . 
0<k<N I ^— ' \ \/ 5n It 



\ 



2p 



yields that F is Lipschitz continuous so using Lemma [3] we show that Iq < j£zp- Next, we have that 



h < C E 



max } / h(Y s )ds - S N h(Y t 



3=0 " L i 



J 2p 



4 p e 



k-l 



max V/i F t . - h(Y t N ) 

3=0 



2p 



On one hand, thanks to assumption (7^|T|) and Lemma [31 



4 p e 



k-l 



^ a /JE^)-M? t f) 

3=0 



2p 



N-l 



< CS N E 

3=0 



0<k<N 

On the other hand, using an integration by parts formula, 



\h(Y tj )-h(Y») 



2p 



< 



c 



h := E 



= E 



/ h(Y s ) - h(Y tj )di 



max 

0<k<N 



3=0 " l i 
k-l r t j+1 



53/ (t J+1 -«)((&&' 
,=o •'ti v 



er 2 /l" \ 2 P 

+ —)(n)d« + ^ / (y.)tw 4 



< 2 2p - x E 



max 

0<fc<A r 



tfc cr 2 h" 
(t s - s)(bti + —)(Y s )ds 



2p- 


+ E 








max 


L 






0<k<N 





(r s - s)cr/i'(r s )dW a 



2p 
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where we denoted by t s the lowest discretization point greater than s : t s = [j 2 -]^. Using Jensen's 
inequality for the first integral and the Burkholder-Davis-Gundy inequality for the second, we obtain 



h < C E 



max t 

0<k<N 



2p-l 



(t s - s) 



2p 



{bh' + —){Y s ) 



2p 



ds 



(T s -s) 2 \ah'(Y s )\ 2 ds 



< 



C 
N 2 p 



E 



2p 



\ah'(Y s )\ 



2p 



ds. 



Under the assumptions of Lemma [31 sup 0<t<T E{\Y S \ 2p ) < oo (see Problem 3.15 p. 306 of Karatzas and 
Shreve [1991] for example) so, with the help of assumption (?42]), we conclude that Ii < and hence 
I\ < jfz? ■ We now turn to the last term. Using Burkholder-Davis-Gundy inequality, we get 



h < CS P N E 



In -i 



\j=o 



HY s )ds- 



A 



>(Y. N ) rtj+i 



>N 



(w s -w tj )ds vv> 



N-l 

< S N E 



\ 



>N 



(W a -W tj )d8 VV 



2>r 



(8) 



Assumption yields that the two terms appearing in the square root are bounded from below by ip > 
so we have that 



N-l 
j=0 



^J t h+1 m)ds- U(Y t f) + 



aiP'(Y t f) r* i+ i 



?N J t 





2p" 







N-l 
N-l 

3=0 



^(Y s )ds - \^>(Y»)8 N + (n//(Y») I {W s - W tj )d 



2;/ 



where 



and 



ll = E 



iP(Y s )ds - ip(Y tj )S N + o^'(Y t] ) / (W s - W t] )ds 



2p 



J|=E 



S N U(Y tj ) - ^{Y»)\ + Unl/QTt,) - <nf/(Yg)) / (W s - W tj )d 



2p 



Again, integrating by parts yields that 



(9) 



1 2 = E 



(t j+1 - s) (ai>'(Y s ) - <jyj'(Y tj ))dW s + {(in// + —^'){Y B ))di 



2p 
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We control the stochastic integral term as follows 



(tj+i - s)W(Y s ) - ai>' \Y tj ))dW s 



2p 





rtj+i 










<C5 P N ^E 


/ (*j 

Jtj 


f i - *) 2p 




a^'{Y tj )\ 2 Us 


_ 

<CS%^ j 


ftj + l r 

' E 






\ 2P ] ^ 




ftj + l r 

' E 

tj 1 






ds 






/■tj+i 
' 










<C5% 













The third inequality is due to assumption (7-43J and the fourth one is a standard result on the control of the 
moments of the increments of the solution of a SDE with Lipschitz continuous coefficients (see Problem 3.15 
p. 306 of Karatzas and Shreve [1991] for example). 

We also control the other term thanks to assumption (T-^ : 



(t j+1 -s)(b^' + —i>")(Y s )ds 



2// 



< 




f 












"tj+i 


< 






3 


< 


C6% 





{t j+ i-B)*»\w + -r)(y.)?*d8 



2p 



ds 



Hence, 1 2 < T^xp- To conclude the proof of the theorem, it remains to show a similar result for 1^ 



< C[S 2 N p E 



Sn (^)- mf )) + (^'(Y tj )- ^'(Y t f)j / (W s - W t] )d, 



2p 



N 



Y t . - Yf. 



2p 



3p 



Yu ~ Y t 



N 



2p 



< 



c 

N 4 p 



The second inequality is due to the fact that ip is Lipschitz continuous (thanks to assumption (7411)) f° r the 
first term and to the independence of (cnp'iYtj) -aip'(Y t f) \ and J^ +1 (W s - W tj )ds for the second term. □ 



Remark 5 — Our scheme exhibits the same convergence properties as the Cruzeiro et al. [2004] scheme. 
Apart from the fact that it involves less terms, it presents the advantage of improving the multilevel Monte 
Carlo convergence. This method, which is a generalization of the statistical Romberg extrapolation method 
of Kebaier [2005], was introduced by Giles [2008b, 2008a]. 

Indeed, consider the discretization scheme with time step 5 2 n — to ■' 



V0 < k < 2N- 1, X 



2N 



= X™ + p( F{Yl N +1)T ) F(Y™)) + 5 2N h(YM) + VT 



A 



4>(Y™) + 

2N 



^'(YM) 

2N 

$2N 



(fc + l)T \ 

/ 2N (W s - Wkr)ds J VV> (Bi 



B kT ) 

2JV / 
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Denote by vi N = yl — p 1 



\ 



2N\ 

kT } 
IN 



JkT 2N (W a — W kT )ds V ip the random variable which 



multiplies the increment of the Brownian motion \B{k+i)T — Bwr_ ) . Because of the independence properties, 

V 2n — 2 AT / 

~ \ (~ N \ ~ N 

X t N ) has the same distribution law as the vector ( X+, } defined inductively by X tQ = log(so) 



and 



0<k<N 



0<k<N 



■N 



VO < k < N - 1, X tk+i =X tk+P IF(Y») - F{Y»)\ + 5 N h{Y«) 



whe 



Ai?f +1 = V2 



2N 
2k 



\ 



>N 



B, 



B lkT 
2N 



V, 



2N 
2k+l 



B, 



B, 



>4n 



l! 2N ) 
U 2k+1! 



(10) 



Going over the proof of the theorem, one can show in the same way that 



0{N~ 2 ) 





~N 


2" 


max 


x t -x? N 




0<k<N 













(11) 



Hence, one can apply the multilevel Monte Carlo method to compute the expectation of a Lipschitz continuous 
functional of X and reduce the computational cost to achieve a desired root-mean- square error of e > to a 
0(e- 2 ). 

As a matter of fact, the particular structure of our scheme enabled us to reconstruct the coupling which 
allows to efficiently control the error between the scheme with time step jf and the one with time step -£r. 
This does not seem possible with the Cruzeiro et al. [2004] scheme. 

From a practical point of view, it is more interesting to obtain a convergence result for the stock price. 
It is also more challenging because the exponential function is not globally Lipschitz continuous. We can 
nevertheless state the following corollary with some general assumptions and we will see in the next section 
that we can make them more precise in case (^t)«e[o,T] 1S an Ornstein-Uhlenbeck process. 



Corollary 6 — Let p > 1. Under the assumptions of Theorem^ and if 
(HI) 



3e > such that E 



E 



0<k<N k 

then there exists a positive constant C independent of N such that 



max e 

0<k<N 



(2p+e)X* 



< oo 



E 







2p 


max 


e Xt k - e ** 




0<k<N 







< 



c 
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Proof : Using Holder inequality we have that 



E 


max 




2p" 


< E 


max 




0<fc<Af 








0<fc<JV 



2p 



< E 



max 

0<k<N k 



max e 

0<k<N 



2p 
2p + c 



x E 



max 

0<fc<JV 



N 



We conclude by assumption (?4Zl) and Theorem [2] 



Remark 7 — _ffad we introduced a new cut-off to our scheme as follows 



\ 



assumption (h^Tty would have been induced by assuming that the functions F, f and h are bounded. 



1.2 Special case of an Ornstein-Uhlenbeck process driving the volatility 

For many stochastic volatility models, the process (^t)tg[o,T] which drives the volatility is an Ornstein- 
Uhlenbeck process. For example, this is the case for all the models cited in the introduction but the Heston 
model. Therefore, it is useful to focus on this particular case. We will hereafter suppose that (^)tg[o,T] is 
solution of the following SDE 

dY t = vdW t + k(9 - Y t )dt (12) 

with v, k and three positive constants. Since exact simulation is possible, we can replace the Milstein 
discretization by the true solution in our previous scheme : 



WeakTraj_l scheme when Y is an O-U process 



K + i = X Z + P (F(Yt k+1 ) - F(Y tk )) + 6 N h(Y tk ) 



>N 



(W s -W tk )ds AB fc+ i 



(13) 



Note that we require the exact simulation of both (Yt fc , Yt k+1 ) and f t k+1 (W s — Wt k )ds. The unique solution 
of CE2) is F t = y e- Kt + 0(1 - e~ Kt ) + v f* e" K( *" s) dTU s and one can easily deduce that, Vfc e {0, . . . , N - 1}, 



Y tk+1 -e-^Y tk 
J^(W s -W tk )ds 



• Af (M, T) 
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/ 8(1 - e-riv) \ / |_(i- e -2«5«) %(l-e-" s »(l + KS N )) 

where M = I and T = 

V J \ £(1_ 6-^(1 + ^)) & 

We first state the following technical lemma whose proof is postponed to the appendix : 



Lemma 8 — V c x > 0, c 2 E [0, 1) 



E 



^ e ci sup 



< t <T|V t | 1 + C2 



< oo. 



As might be expected, it is possible to weaken the assumptions of Theorem O In particular, we relax 
the assumption on the lower bound of the volatility (742) an d replace it with a weaker one (see assumption 
CHM below). 



Theorem 9 — Let p > 1. Suppose that Y is solution of (ltfy and that the scheme is defined by (J 
assumption (T~^) of Theorem [H and if 

(H8) f is a C 3 function 

(H9) there exist three constants cq > 0,ci > and c 2 £ [0, 1) such that, Vy £E R, 

<cne c ^ 1+c2 



Under 



K (9-y)h'{y) + -h"{y) 



\h'(y)\ <c e c ^ 



l + c 2 



ip"(y)\<c e c ^ 1+C2 



< c e ClMl 



then there exists a constant C independent of the number of time steps N such that 



max 

0<k<N 



N 



2// 



< 



c 



The same result holds true when we replace assumption (T~t^) by 
(HIO) There exist two positive constants C and e such that Vy € R, 

<P(y) > o 

supE(V (1+e) (^)) < oo 

t.<T V / 

< OO. 



t<T 

sup E 

t<T 



Proof : The proof of the first part of the theorem repeats the proof of Theorem [2] with fewer terms to 
control because of the exact simulation of (Yt)te[o,T]- At the places where we used assumptions (7421) and 
(74U), we use assumption (Tf^ together with Lemma [51 
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We now focus on the second part of the theorem. According to equation (JSJ), all we have to show is the 
existence of a positive constant C independent of N such that Vj € {0, . . . , N — 1} 



7- r +i - 



S N Jt 3 



\ 



ON 



(W a -W tj )ds y± 



2p- 



< 



c 

N 2 p 



We will adopt the following notations 

- Dj = ^(Y t] ) + ^1 (W s - W tj )d»)v± 

Thanks to assumption ff JTU|) , we have that Vj € {0, . . . , N— 1}, Aj > and Dj > 0. The idea of the proof is 
to isolate the case where Dj is small which is problematic since the square root is not Lipschitz continuous 
in the neighborhood of : 



'A 3 - VDj 



2p 



'Aj - y/Di 



2p 



1 



{D j <i,(Y t] )/2} 



'Aj - ^Dj 



2l> 



1 



{D j >i,(Y ti )/2} 



2P 



+2 2p- 



MD J <4,(Y t] )/2} 
2 p 



A p j ^(Y tj ) 



A j - D j\ 2Pl {D j >^(Y tj )/2} 



We take the expectation and apply Holder inequality to obtain 

2pl 



E 



'Aj - ^Dj 



with 



and 



ei = E 



62= E 



2p 



l+e 



2" 



A v - r(Y tj ) 



Let us begin with the second term. Following the estimation of l\ m the proof of Theorem [21 we show that 



Aj-Dj\^ 



= E 



< CS 



3 AT Jt 



— TP(Y a )ds- [ij(Y tj ) + 





2pl±i- 







2p^ 
N 



Thanks to assumption ("HUB and Jensen's inequality, we also have that 




2'-' 



l+e 



^ p {Y tj ) 




< C 



ll>p(.l+e)(Y s ) 



ds + 2 p ( 1+£ )E 
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Hence ti < jpp- Now let us turn to t\. Note first that assumption ff/fTU]) enables us to show that there 



exists a positive constant C independent of TV such that ( E 

is left to prove is that P < — -p 2 -^ < 
Va > 0, 3C a > such that P (Dj < 



< C. Finally, what 



N 



In fact, we can show that 



< 



N° 



P 



uib'(Y t .) f^ +1 , s ib(Y t .) 

^ V *> } ' {W s -W tj )ds<-^^- 



1N 



= P I IGI > 

where G is a centered reduced Gaussian random variable independent of Y tj 

V3i>(Y t 



Thanks to assumption (Ti HO]) . 3C > s.t. P MG| > 2 ^l^]^ t y\ < 2P ( G > and usin S the 

following standard upper bound of the Gaussian tail probability : Vi > 0, P(G > t) < - j= , we conclude. □ 



Remark 10 - 

• The fact that we can simulate exactly the volatility process without affecting the order of convergence 
of the scheme is yet another advantage of our approach over the Cruzeiro et al. [2004] scheme. On 
the other hand, the Kahl and Jackel [2006] scheme allows the exact simulation of (^t)tG[o,T]- Applied 
to the SDE |3J), it writes as 

vUK _ viJK , ( _ f 2 ( Y t k+1 ) + f 2 (Yt k ) \ . \ AW 

X t k+1 ~ X t k +{r 1 J °N + Pf{Yt k )/\W k+ i 

+ V 1 - P ^ AB k+i + y/ [Ytj [(AWt+i) -5n) 

Note that it is close to our scheme insofar as it takes advantage of the structure of the SDE (for example, 
unlike the Cruzeiro et al. [2004] scheme, it allows the use of the coupling introduced in Remark^. The 
main difference, which explains why our scheme has better weak trajectorial convergence order, is that 
we discretize more accurately the integral of f(Y t ) with respect to the Brownian motion (B t )t£[o,T] ■ If 
instead of a trapezoidal method, one uses the same discretization as for the WeakTraj-1 scheme, then 
it can be shown that this modified UK scheme will exhibit a first order weak trajectorial convergence. 

• One can easily check that this theorem applies for the Scott's [1987] model (and therefore for the Hull 

2 2y 

and White [1987] model) where we have h{y) = r — ^| pa o e y (^(0 — y) + |) and ipiu) = °~o e2v ■ 

The Stein and Stein [1991] and the quadratic Gaussian models do not satisfy the assumption \ip'{y)\ < 
Cv{y). 

• It is possible to improve the convergence at fixed times up to the order I . Following Lapeyre and Temam 
[2001] who approximate an integral of the form f t h+1 g(Y s )ds for a twice differ entiable function g by 

S N g(Y tk ) + vg'{Y th ) f£ +1 (W s - W tk )ds + (k(9 - Y tk )g'{Y tk ) + £g»(Y th )) S -f, we obtain the following 
scheme 

OU -Improved scheme 

X t N k+1 =X t N k +P (F(Y tk+1 ) - F(Y t J) +h k + yr^V^fe AB k+1 (15) 



14 



where h k = 5 N h(Y tk ) + vti{Y tk ) J^ +1 {W S - W tk )ds + (k(6 - Y tk )h'{Y tk ) + ' (Y tk )) A -f and 

!^lj^(W s - W tk )ds + (k(6 - Y tk W{Y tk ) + ^"{Y tk )Y-f) V ±. 
Mimicking the proof of Theorem fJl one can show that 

2" 



max E 

0<k<N 



= O (N- 3 ) 



where Xt k and Xfr have respectively the same distribution as Xt k and X t 



N 



X tk =X + p{F{Y tk ) - F(y )) + J " h{Y s )ds + ^fl- P 2 J ~ j " ^(Y s )ds B t 



and 



fe-i 



X" = X Q + P (F(Y tk ) - F(yo)) + V 7 !^? 



fe-i 



\ 3=0 



As for the stock, we can prove the same convergence result under some additional assumptions which 
are more explicit than assumption of Corollary [5J To do so, let us make the following changes in our 

scheme so that we can control its exponential moments : 



X t N k+1 =X t N k+P (F(Y tk+1 ) ~ F(Y tk )) + S N h(Y tk ) 



Sn 



(16) 



(W s -W tk )ds) A ip(Y tk ) V tp AB 



k+l 



Proposition 11 — Suppose that Y is solution of ilty) and that the scheme is defined by hi 6(1 . 
Under the assumptions (H\Bj), (H\§j) and (HI 0(1 of Theorem^ and if 

(Ull) there exists (3 € (0, 1) and K > such that Vy £ R 

\h(y)\ + \F(y)\ + \f'(y)\ < K(l + \y\^) 

\f(y)\ < K(l + \yf) 

then, Vp > 1, there exists a positive constant C independent of N such that 



E 



max 

0<k<N 



e x <> 



x 



2p 



< 



c 



The same result holds true if one replaces assumption (Hi 0(1 by assumption (HJ^) together with the assumption 
that 3C > for which Vy € R, \i>'(y)\ < Cip(y). 



Proof : We go over the proof of Corollary [6] The fact that E 



max < fc <7v 



X, 



XI 



N 



Ap 



is not a straightforward consequence of Theorem [5] anymore because we have introduced some changes in 
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our scheme. However, looking through the proof of the theorem, one can see that it is enough to prove the 
following inequality : Vj G {0, . . . , N — 1} 



Sn J i 



^(Y s )ds 



\ 



{W s -W t] )ds Ai>(Y t .)V± 



2,f 



< 



c 

N 2 p 



(17) 



When ip is hnite, since jj_ J t 3+1 i()(Y s )ds is smaller than ^(Y tk ) = tp, we can remove the new cut- 
off from the left hand side of (IT7|) and then proceed like in Theorem [HI When t/j = +oo, on the event 
(^ip{Y tj ) H — - ' g N tj - J* t * 3+1 (W s — Wtj )ds*j < ^(Ytj), we recover our original scheme and we prove (|17l) like in 

Theorem [9] Then, using the Gaussian arguments developed in the end of the proof of Theorem [9l we control 
the probability of the complementary event to conclude. 

Now, what is left to prove is that assumption (?47]) is satisfied. On the one hand, we have that 

4p~ 



max S f p 

0<k<N k 



max Sn 

0<k<N 



n S*ds 



f{Y s )S s (pdW s + y/1 - p 2 dB^ 



C I 



(Sf p (l + f 4p (Y t )) 



dt 



< C 1 + 



Thanks to assumption C hHU\f and LemmafHl there exists C > such that ^/E((l + f ip (Y t )) 2 ) < C. Observe 
that conditionally on (Yt)te[o,T]i 



X t ~Af[ log(so) + p(F(Y t ) - F(y )) + / h(Y s )ds , (1 - p 2 ) / f(Y s )ds 



(18) 



so, by Jensen's inequality and assumption ff fllll) 



= E ^ e 8p(log(^o)+p(i ;, (y t )--F(?/o))+/ t '»(n)^) e 32p 2 (l-p 2 )/ t / 2 (Y s )d s N 

< £ ^ e 8rtlog(*o)+p(^(i' t )-F(ao)))Iyr t e <8p/ l (y s )+32p 2 (l-p 2 )/ 2 (n)) ds ^ 

< CE ^ e Csu Po<t<T l^tl 1 " 



Using Lemma HI we deduce that E 



ma,x 0<k<N S t p 



< oo. 



On the other hand, using Cauchy-Schwartz inequality, we have that 



max e *» 

0<fc<JV 



fe-i k—i 
m«c exp I Ap I X + p(F(Y tk ) - F(y j) + ]T 6 N h(Y tj ) + £ Jl^f 

3=0 j=0 



\ 



Wit 



(W s - W tj )ds A V(y tj ) V ip ABj + i 



where 



E" = E 



max e 

0<fc<JV 



Sp(X +p{F(Y tk )~F(y a ))+J2 k jZo s Nh(Y tj )) 
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and 



E» = E 



max e 

0<k<N 



< oo. 



Using the same argument as before, we show that E^ < CE fe Csup °^ T \ Yt \ 1+f> 

Denote by Dj = [^{Y t] ) + ^J' 3 f t 3+1 {W s - W tj )ds) A ^)V^. Using Doob's maximal inequality 
for the positive submartingale ^e 4p ^~ 



-p 2 E j=0 V D i AB i+i / see Theorem 3.8 p. 13 of Karatzas and 

/ 0<k<N 

Shreve [1991] for example), we also have that 



E? < 4E 



-p 2 E"=o \M ab j+i 



'N-l 



4e n 



,32 P 2 S N (l-p 2 )D;j 



j=0 



< 4E I max e J 

i 0<fc<iV-l 



By virtue of assumption f? dll[) , E% < 00 which concludes the proof. 



2 A second order weak scheme 

Integrating the first stochastic differential equation in (0} gives 

X t = log(flo) + p(F(Y t ) - F{y a )) + f h(Y s )ds + yJX-p 1 f f(Y s )dB s (19) 

Jo Jo 

We are only left with an integral with respect to time which can be handled by the use of a trape- 
zoidal scheme and a stochastic integral where the integrand is independent of the Brownian motion. Hence, 
conditionally on O^tgro.T]) 

X T ~N (log(so) + P(F(Y T ) - F(yo)) + m T , (1 - p 2 )v T ) (20) 

where vtlt = Jq h(Y s )ds and vt — J n T f 2 (Y s )ds. This suggests that, in order to properly approximate the 
law of Xt, one should accurately approximate the law of Yr and carefully handle integrals with respect to 
time of functions of the process (Yt)t<z[o,T\- We thus define our weak scheme as follows 

Weak_2 scheme 

X? = log( So ) + P (F(Yt) - F(y ))+m^ + ^{l-p 2 )v^G (21) 

, — N x ^N-l ^fj+M^f ) -N X f 2 ^ k )+f{y^ k+1 ) lTP N . 

where mif. = o N )^ k=0 - — 5 = o N }^ k=0 - — 5 !tt H K Y t k )o<k<N is the Nmomiya- 

Victoir scheme of (Y t ) t ^[o,T] and G is an independent centered reduced Gaussian random variable. Note 

TV JV . . . N 

that, conditionally on (Y tk )o<k<N , X t is also a Gaussian random variable with mean log(so) + p(F(Y T ) — 
F(y j) +Wit and variance (1 — p 2 )vt ■ 
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It is well known that the Ninomiya and Victoir [2008] scheme is of weak order two. For the sake of 
completeness, we give its definition in our setting : 

\ VO < k < N - 1, Yl +i = exp exp ((W tk+1 - W tk )V) exp (£?V Q ) (<) 

where Vq : x i— > — ^aa'(x) and : a; i— > cr(x). The notation exp(tV)(x) stands for the solution, at time t 
and starting from x, of the ODE rj'{t) = V(j](t)). What is nice with our setting is that we are in dimension 
one and thus such ODEs can be solved explicitly. Indeed, if £ is a primitive of y : ((t) = J Q y^ds, then 
the solution writes as r](t) = (t + C{x)). 

Note that our scheme can be seen as a splitting scheme for the SDE satisfied by 
{Z t = X t -pF{Y t ),Y t ) : 



dZ t = h(Y t )dt+^T~p^f(Y t )dB t 
dY t = b(Y t )dt + a{Y t )dW t 



(22) 



The differential operator associated to (|22|) writes as 

Cv(z,y) = h{y)- + b(y)- + -f 1 ^ + \v)q^ = C Y v(z,y) + C z v(z,y) 

where C Y v(z,y) = b(y)$g + ^0 and C z v(z,y) = h(y)§£ + ^lf(y)^. One can check that our 
scheme amounts to first integrate exactly Cz over a half time step then apply the Ninomiya- Victoir scheme 
to Cy over a time step and finally integrate exactly Cz over a half time step. According to results on 
splitting (see Alfonsi [2009] or Tanaka and Kohatsu-Higa [2009] for example) one expects this scheme to 
exhibit second order weak convergence. We will not use this point of view to prove our convergence result 
stated in the next theorem, since we need to apply test functions with exponential growth to Xt to be able 
to analyse weak convergence of the stock price. 



Theorem 12 — Suppose that p G (—1, 1). // the following assumptions hold 

(HI 2) b and a are respectively C 4 and C 5 , with bounded derivatives of any order greater or equal to 1. 
(HI 3) h and f are C and F is C 6 . The three functions are bounded together with all their derivatives. 
(H14) ±>0 

then, for any measurable function g verifying 3c > 0,/i G [0,2) such that Vx G R, |.g(a;)| < ce' 2 ^, there 
exists C > such that 



(g{X T ))-E(g{X?j) 



C_ 

~ N 2 



In terms of the asset price, we easily deduce the following corollary : 



Corollary 13 — Under the assumptions of Theorem \1 6 A for any measurable function a verifying 3c > 
0,/i G [0,2) such thatyy > 0, \a(y)\ < cel log ^l", there exists C > such that 



E(a(S T )) - E I a(e 



< 



C_ 

iV2 
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Proof of the theorem : The idea of the proof consists in conditioning by the Brownian motion which 
drives the volatility process and then applying the weak error analysis of Talay and Tubaro [1990]. 

As stated above, conditionally on (Wt)te[o,T\, both Xt and X T are Gaussian random variables and one 
can easily show that 



E 



g{X T )-g(x") 



g(x)E 



exp 



( (x-\og(s a )+pF(y ] 



)-pF(Y T )-m, T y 

-P 2 )V T 



exp 



(x-\og(s )+pF(yo)-pF(Y T )-m%)< 



2tt(1 - p 2 )v% 



dx 



For x £ R, denote by j x the function 

~f x : R x R x R* -> 

(y, m, v) i-> 



(x-\og(s„)+pF(y a )-pF(y)- m y 
eX P I 2(l-pi)u 



y/2n(l - p 2 )v 



so that e < f R g(x) E ^y x (Y T , m T , v T ) - 1x(Yt > ™t > w t ) 
following intermediate result : 3C, K > and pg N such that 

V.t € R, |E [ lx (Y T ,m T ,v T )- lx (Y^,m^,v^)\ < ^e"**" "( I |.r "). 

We naturally consider the following 3-dimensional degenerate SDE: 

f dy t = <7(Y 4 )dW t + 6(y t )dt; Y = y 

dmt — h(Yt)dt; mo = 

dih = f 2 (Yt)dt; v = Q 

Note that (Y^ ,rn^ ,v^) is close to the terminal value of the Ninomiya-Victoir scheme applied to this 3- 
dimcnsional SDE. In order to prove (I23p . we need to analyse the dependence of the error on x and not only 
on N. That is why we resume the error analysis of Ninomiya and Victoir [2008] in a more detailed fashion. 
For x e R, let us define the function u x : [0, T] X R X R X R^ -> R by 



dx. Consequently, it is enough to show the 
C '• - (23) 



(24) 



u x (t,y,m,v) 



E 



x {{Y 



T-t,m T -t,VT-t 



\(y,"i,v] 



where we denote by (Ir-t, wir-t 5 VT-tY y ' m ' v ^ the solution at time T — t of (|24|) starting from (y, m, v). 

The remainder of the proof leans on the following lemmas. We will use the standard notation for 
partial derivatives: for a multi-index a — (a\, . . . , ay) € N d , d being a positive integer, we denote by 
|a| = ax H h ay its length and by d a the differential operator 9'"' /d" 1 . . . d^ d . 

Lemma 14 — Under assumptions (T- H12)) . (T~ M3\) and (h\lj$, we have that 

i) u x is C 3 with respect to the time variable andC 6 with respect to the space variable. Moreover, it solves 
the following PDE 

{ d t u x + Cu x = 

{ (25) 
[ Ux{T,y,m,v) = j x (y,m,v) 

where C is the differential operator associated to 
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ii) For any multi-index a £ N 3 and integer I such that 21 + \a\ < 6, there exists C^ a ,Ki a > and 
(j>i,a,Qi,a) G ^ 2 such that 

V(t, y, m, w) G [0, T] x D t , |d|d a u x (t, y, m, v)\ < C l , a ?r K *»'* (1 + |a;|*'.«) (1 + |y|*.» ) 

where T> t is the set R x [— isup zgR |/i(z)|,isup z6IR \h(z)\] x [tip,tip\. Note that ip and ip are finite by 
virtue of assumptions f H13\) and (1- \1J$ . 



Lemma 15 — Under assumption (Hlty), 



Vg G N, sup E ( Y t 

0<k<N K ' 



< oo 



Now, following the error analysis of Talay and Tubaro [1990], we write that 



N-l 



Jx(Y T , rn T , v T ) - 7 X (Y T , m^^r) < r) k (x) 



k=0 



where n k (x) = E u x (t k +i,Y" ,mf ) - u x (t k ,Y? ,m?,vF) 



V x v^'— 1 

'"t fe _ ° N 2^j=0 2 " iivL "t k — ^j=0 T 

first term in the expectation and Taylor's formula together with PDE (j25|) for the second, we get 



and v tk = 5 N J2j=o 



fc -l / 2 (^)+/ 2 (7g +1 ) 



and V0 < fc < TV, 



Using the Markov property for the 



r]k{x) 



(t>x(tk+i,Y th ,m? k ,v%) - u x (t k+1 ,Y tk ,m^ k ,v^ k ) - 6 N Cu x (t k+1 ,Y tk ,m^ k ,v^ k ) 



— iV 



1 f tk + 1 3 U X . —N _ N _ 



dt 3 



(t,Yt k ,rnl,vl){t~t k ) 2 dt 



where 



<j>x(t k+ i,y,m,v) = E 



M x (t fc+ i,r ti ,m + djv ,w + 0Af 



f(K' v ) + f 2 (y). 



Denote by T y the function z n- u x {t k+ \,z, m + Sn Mfl+Mg) ; w _|_ ^Lifl+LM). Using Taylor's formula we 
can show that Vz G R, 



r„(z) = r Vtl (z) + 5 N v Vt2 {z) + -fr y , 3 ( z ) + r ( z ) 



where 



T yA (z) = u x (tk+i,z,m,v) 

h(z) + h(y) du x f 2 (z) + f 2 (y) du Xrj 

r vA z ) = o -5zr(*fc+i' z ' m ' w J H o sr(*k+i»^ m >«) 



9m 



^) + %)V3 : 



2 «2„ 



-(tfe + i,z,m,w) + 



+2 



2 y ,9m 2 

Mz) + %)/ 2 (z) + / 2 Q/) gV, 
2 2 <9m<9v 



<9v 

/ 2 (z) + / 2 ( y )\ 2 av 



dv 2 



(tk+i,z,m,v) 



(t k+1 ,z,m,v) 
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and 



+ 



\ 2 / h(z)+h(y) \ d 3 ^ (. , m , . hj^+h^y} f f\z)+f 2 (y) 

J \ 2 j 3m9i) 2 ^Hb z ! m + ' 2 , w + r 2 



I q ( f 2 {z)+f A {y) Y ( h(z)+h(y) \ _& 

~r° \ 2 ] \ 2 ) dmdv 2 



So, 



E 


r„,i 







<5jvE 



rj,,2(V'ti' ! 'y 



— E 



</>x,l(tfc + i,t/,m,f) 0i,2(tfc+i,3/,m,w) 



(26) 



(27) 



With a slight abuse of notations, we define the first order differential operators Vb and V acting on C 1 
functions by Vq£(x) = Vo(x)£'(x) and F£(x) = V(x)£ > '(x) for £ e C 1 (R). We make the same expansions as 
in Ninomiya and Victoir [2008] but with making the remainder terms explicit in order to check if they have 
the good behavior with respect to x. We can show after tedious but simple computations that 



c,i{tk+\,y,m,v) = 



+ (4V 2 r\,,i(y) + W V 2 r yA (y) + 2V 2 V r yA (y) + V 4 T yA (y)) + E (J2i(i/)) 



v,3(t k+ i,y,m,v) 



S N T y , 2 (y) + -f (V 2 T ya (y) + 2V T ya (y)) + E (R 2 (y)) 



T y , 3 (y) + E(R 3 (y)) 
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where 

Ri(y) 



/</ So" Jo 2 V 3 r y Ae S3Vo e w ^ v e^ v o {y))ds3ds2dsi 

+ C 6N Io 1 Jo* Jo* Io* Jo" V e T y Ae^ v e^(y))d S(i d S5 d S4 d S3 d S2 d Sl 

+ ¥C N Jo" Jo 2 Io' V*V r y Ae s * v e^o {y))ds4ds3ds2dsi 

+ fC SN Jo" V*Vo 2 r yil (e s > v e^(y))d S2 d Sl 

+ Io* Io' Io 2 V 3 T y Ae S3V °(y))ds 3 d S2 d Sl + 5 -f // / Sl V 2 V 2 T yA (e s ^> (y))ds 2 d Sl 

+f j-T ^% 4 ( e ^(y))^i + 4*/^ / o Sl Vb 3 r y a(^ y °(2/))^i 



i? 2 (y) = M/ 3 / sl Vb 2 r a ,2(e S2V °e^^^(y))d S2< i Sl 

+ C™ Jo" IT Jo 3 V^ 2 (e^e^>(y))d S4 ds 3 d S2 d Sl 

+ S itC 5N I" V 2 V Q Y ya {e^e^ {y))ds2dsi+ J^- £ V 2 T y Me^(y))ds 2 d Sl 
+ % L io" T " V V 2 T y . 2 {e^( y ))d Sl + 5 -f VQ 2 r y)2 (e s ^"(y))d Sl 



+ J Q 2 v r y ^ v °( y ))d Sl 

Putting all the terms together, one can check that 



S 2 

<j) x (t k+1 ,y,m 7 v) = u x (tk + i,y,m,v) + S N £u x (t k +i, y, m, v) + -^-C u x {tk+\, y, m, v) + R(y) 



where R(y) = E 



R (Yf i ' v )+R 1 (y)+R 2 (y)+R 3 (y) 



N-l 



E 



j x (Y T ,mT,v T ) -Tx^T.my,^) < ^ E 



fe=0 



Finally, 

1 f tk+1 d 3 u x 



{t,Y tk ,mZ,vZ){t-t k ydt 



(28) 



■A 1 , 



R(Yt, 



From Lemmas [14] and [T5l we deduce that there exists C± , K\ > and p\ G N such that 

' ' r 1 d 3 U, r , „. 



fc=0 



at 3 



■(t,Y't k ,mf k ,v^ k )(t-t k ) 2 dt 



<^c ie - K ^ 2 (i+\xn 



(29) 



is of order 



On the other hand, a close look to ([26]) and ([28]) convinces us that the term E 

and that it involves only derivatives of u x and of the coefficients of the SDE ([24]) . So, thanks Lemmas [T4l 
and[l5l there exists C 2 , K 2 > and p 2 £ N such that 

JV-l 



fc=0 



^^e-^ 2 (i + ixr) 



(30) 



From (l29l) and (f30l) we deduce the desired result ([23l) to conclude 



□ 
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Remark 16 — 

• The theorem does not cover the case of perfectly correlated or uncorrelated stock and volatility which is 
not very interesting from a practical point of view. 

• As for plain vanilla options pricing, observe that, by the Romano and Touzi [1997] formula, 

E(e- rT a(S T )\(Y t ) m )=BS a , T L e KF(r T )-F( 2/0 )) +mT+ (<i^-,)T ) 

where BS at T(s, v) stands for the price of a European option with pay-off a and maturity T in the Black 
& Scholes model with initial stock price s, volatility y^v and constant interest rate r. When, like for 
a call or a put option, BS ay r is available in a closed form, one should approximate E (e~ rT a(ST)) by 

1 M 

where M is the total number of Monte Carlo samples and the index i refers to independent draws. 

Indeed, the conditioning provides a variance reduction. We also note that what is most important is to 
have a scheme with a high order weak convergence on the triplet (Y t ,m t ,v t ) te [ T -\ solution of the SDE 
\2$ , which is the case for our scheme. 

• In the special case of an Ornstein- Uhlenbeck process driving the volatility (i.e {Yt)t£[o,T} is solution 
of the SDE \12]) ). one should replace the Ninomiya-Victoir scheme by the true solution. We can 
then prove more easily the same weak convergence result: at step of the preceding proof, we 
apply ltd 's formula instead of carrying out the Ninomiya- Victoir expansion. Moreover, we can prove, 
following the same error analysis, that the OXJ -Improved scheme 115\) also exhibits a second order weak 
convergence property. Better still, it achieves a weak trajectorial convergence of order | on the triplet 
(Y t , m t , v t )t£[o,T] which allows for a significant improvement of the multilevel Monte Carlo method, as 
we shall check numerically. 



s e 



p(F(Yt 



)-ffeo))+m"' i + ( 



r)T ~ P 2 ) V 1 



3 Numerical results 

For numerical computations, we are going to consider Scott's model We use the same set of param- 
eters as in Kahl and Jackel [2006] : So = 100, r = 0.05, T = 1, <t = 0.25, y = 0, k = 1, 6 = 0, v = = 
-0.2 and / : y h> cr e y . 

We are going to compare our schemes (WeakTraj_l, Weak_2 and OUJmproved) to the Euler scheme with 
exact simulation of the volatility (hereafter denoted Euler), the Kahl and Jackel [2006] scheme (UK) and 
the Cruzeiro et al. [2004] scheme (CMT). 



3.1 Numerical illustration of strong convergence properties 



In order to illustrate the strong convergence rate of a discretization scheme X , we consider the squared 
L 2 -norm of the supremum of the difference between the scheme with time step jj and the one with time 
step 



2N 



E 



max 

0<fc<A' 



X 



N 



x 



2jV 



(31) 



This quantity will exhibit the same asymptotic behavior with respect to N as the squared L 2 -norm of 
;rence b 
Alfonsi [2005]). 



the difference between the scheme with time step S and the limiting process towards which it converges (see 
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In Figure [TJ we draw the logarithm of the Monte Carlo estimation of (f3"Tj) as a function of the logarithm of 
the number of time steps. The number of Monte Carlo samples used is equal to M = 10000 and the number 
of discretization steps is a power of 2 varying from 2 to 256. We also consider the strong convergence of the 

schemes on the asset itself (see Figure [5]) by computing E 



maxo<fe<Ar 



The slopes of the regression lines are reported in Table [TJ We see that, both for the logarithm of the 
asset and for the asset itself, all the schemes exhibit a strong convergence of order \ . Our schemes only have 
a better constant. 
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Figure 1: Strong convergence on the log-asset 



Figure 2: Strong convergence on the asset 
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-0.98 


-0.84 


Asset 


-1.01 


-0.91 


-0.95 


-0.88 


-0.95 


-0.85 



Table 1: Slopes of the regression lines (Strong convergence) 



3.1.1 Weak trajectorial convergence 

Nevertheless, as explained in Remark [5J for the scheme with time step j?, one can replace the increments 
of the Brownian motion (-B t ) t6 [ ,T] by a sequence of Gaussian random variables smartly constructed from 
the scheme with time step This particular coupling is possible whenever the independence structure 
between (B t ) tf =[o t T] an d (Yt)t e [o,T] is preserved by the discretization of the latter process, which is the case for 
all the schemes but the CMT scheme. So we carry out this coupling and we repeat the preceding numerical 
experiment. The results are put together in Figures [3] and U] and in Tabled 

As expected, we see that the WeakTraj_l and the OUJmproved schemes exhibit a first order convergence 
rate whereas the other schemes exhibit a \ order convergence rate. Note that the CMT scheme has a weak 
trajectorial convergence of order one but it is much more difficult to implement the coupling for which the 
convergence order is indeed equal to one. 
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Figure 3: Weak trajectorial convergence on the Figure 4: Weak trajectorial convergence on the 

log-asset (with coupling) asset (with coupling) 
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Table 2: Slopes of the regression lines (Weak trajectorial convergence) 



3.1.2 Convergence at terminal time 

We consider now convergence at terminal time, precisely the squared L 2 -norm of the difference between 

Z and X. 

N ana 2N 

2 



the terminal values of the schemes with time steps and tj?t 



E 



X 



N 



X 



2 A 



(32) 



Note that we introduce a coupling : we write the schemes straight at the terminal time as we did for the 
Weak_2 scheme (see (f2"Tj) ) and we generate the terminal values of the schemes with time steps j£ and X 
using the same single normal random variable to simulate the stochastic integral w.r.t. (-Bt)te[o,T]- Once 
again, it is possible to proceed alike for all the schemes but the CMT scheme. For the latter, we simulate 
the scheme at all the intermediate discretization times to obtain the value at terminal time. 

We also consider the convergence at terminal time of the asset itself. We report the numerical results in 
Figures [5] and [S] and give the slopes of the regression lines in Table [3] 

We observe that, as stated in Remark[T0l the OUJmproved scheme exhibits a convergence rate of order |, 
outperforming all the other schemes. As previously, the WeakTrak_l scheme exhibits a first order convergence 
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Table 3: Slopes of the regression lines (Convergence at terminal time) 
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Figure 5: Convergence at terminal time for the Figure 6: Convergence at terminal time for the 

log-asset asset 

rate. Note also that this new coupling at terminal time improved the convergence rate of the Weak_2 and 
the UK schemes up to order one and, surprisingly, it improved the convergence rate of the Euler scheme up 
to an order strictly greater than the expected ^, approximately 0.67. 

3.2 Standard call pricing 

3.2.1 Numerical illustration of weak convergence 

We compute the price of a call option with strike K = 100 and maturity T = 1. For all the schemes but 
the CMT scheme, we use the conditioning variance reduction technique presented in Remark \W[ 

In Figure [JJ we draw the price as a function of the number of time steps for each scheme and in Figure 
|8]we draw the logarithm of the pricing error : log (|F xact — -P S chcmc|) wnere ^exact ~ 12.82603 is obtained 
by a multilevel Monte Carlo with an accuracy of 5bp, as a function of the logarithm of the number of times 
steps. 

We see that, as expected, the Weak_2 scheme and the OUJmproved scheme exhibit a weak convergence 
of order two and converge much faster than the others. The weak scheme already gives an accurate price 
with only four time steps. The WeakTraj_l scheme has a weak convergence of order one like the Euler and 
the UK scheme, but it has a greater leading error term. Fortunately, its better strong convergence properties 
enable it to catch up with the multilevel Monte Carlo method as we will see hereafter. 

Finally, note that the weak scheme does not require the simulation of additional terms when compared 
to the Euler or the UK schemes. Combined with its second order weak convergence order, this makes the 
Weak_2 scheme very competitive for the pricing of plain vanilla European option. 

3.2.2 Multilevel monte carlo 

Let us now apply the multilevel Monte Carlo method of Giles [2008 b] to compute the Call price. As 
previously, we consider the schemes straight at the terminal time and use a conditioning variance reduction 
technique. We give the CPU time as a function of the root mean square error in Figure [S] (see Giles [20086] 
for details on the heuristic numerical algorithm which is used). 

We observe that both the Weak_2 and the OUJmproved scheme are great time-savers. For the OUJmproved 
scheme, the effect coming from its good strong convergence properties is somewhat offset by the additional 
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Figure 8: Illustration of the convergence rate for 
the call option 



terms it requires to simulate. We can see nevertheless that it is going to overcome the Weak_2 scheme for 
bigger accuracy levels. 




Figure 9: Multilevel Monte Carlo method for a Call option using different schemes 



3.3 Lookback option pricing and multilevel Monte Carlo 

Finally, we consider an example of path-dependent option pricing : the Lookback option. More precisely, 
we compute the price of the option whose pay-off is equal to St — ^oio-te[o,T] St- 
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In order to take full advantage of the good convergence properties of our schemes, we approximate the 
minimum of the scheme by the minimum of a drifted Brownian motion. This is similar to what is done in 
Giles [2008 a]. 

More precisely, for the WeakTraj_l scheme, consider the interval [kjj, (k + 1) j^]- 
Scheme with time step 5 2 n '■ 

We approximate min te [ fe T ^k+i) S] ^* by ™2k ^ ™-2k+i where, V0 < j < 2N — 1, 



I / V 2N / / V 2JV \ ov2N rp 



.2JV 



(l + + f{Y^) ( P (W U+1) ^ - W^) + VT^7(B U+1 ^ - B jA ))) and 
(Uj)o<j<2N-i is an independent sequence of independent random variable uniformly distributed. 

Scheme with time step S jy ■ 

According to Remark [SJ X tk 1 is computed using the Brownian increment AB k+1 given by a linear com- 
bination of iB, 2k+1 \T — B k r_ ,B, k+1 \T_ — B/ a mi t (see (fT0|)). Now, to prevent bias, we are going to 



_ '(2fc+i)^ -^^(k+i)§ --"(2fc+i) w 
approximate min tg [ fe ^ (k+l) T ] &t by the minimum miri[ tfcitfc+1 ] S^" of some Euler scheme S^' ZIV like in the 
scheme with time step To remain consistent, we have to choose 

= A (l + + /(F fc cr ) - WiJ + v^^jj 

In order to ensure a good strong coupling with the scheme with time step <52at, we need to compute 
the intermediate value = e 5 H (l + + ) (p(W tk+1 - W t J + JT=fAB 2k+ ^J us- 

~2JV _ 

ing some Brownian increment AB 2k+l as close as possible to B, 2k+1 yT — B k r_ but such that AB k+1 — 

~2A ~2N ' ~2N 

AB 2k+1 is independent of AB 2k+1 and distributed according to Af(0, 2 tt). Choosing AB 2k+1 of the form 

-B t 1+6 1 B t — B t I / ~2iV 



^=,=5 ^ and maximizing Cov \^AB 2k+1 .B {2k+l) T_ - B k r J = leads 

to a — v 2k + v 2k+1 and b — v 2k+1 — v 2k (see Remark[5]for the definition of v 2N ). 
Finally, we approximate min te j fc r. tj St by rh k A fh k+1 where 



\ - I [ e x ™£ + S e ^ 2N - \ / (e X *% - S e ' 2N \ - 2e 2X *£ f 2 (Y, t ) — 



and 



The numerical results we obtain are very satisfactory. In figure [lOl we draw the CPU time times the 
mean square error against the root mean square error. We see that our schemes perform much better than 
the others. 
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Figure 10: Multilevel Monte Carlo method for a Lookback option using different schemes. 



4 Conclusion 

In this article, we have capitalized on the particular structure of stochastic volatility models to propose 
and discuss two simple and yet competitive discretization schemes. The first one exhibits first order weak 
trajectorial convergence and has the advantage of improving multilevel Monte Carlo methods for the pricing 
of path dependent options. The second one is rather useful for pricing European options since it has a second 
order weak convergence rate. 

We have also focused on the special case of an Ornstein-Uhlenbeck process driving the volatility, which 
encompasses many stochastic volatility models such as the Scott's [1987] model or the quadratic Gaussian 
model. Then, the convergence properties of the previous schemes are preserved when simulating (Y t )o<t<T 
exactly. We have also proposed an improved scheme exhibiting both weak trajectorial convergence of order 
one and weak convergence of order two. 

The numerical experiments show that our schemes are very competitive for the pricing of plain vanilla 
and path-dependent options. Their use with multilevel Monte Carlo gives satisfactory results too. We should 
also mention that the main purpose of our study was the convergence order with respect to the time step. 
It would be of great interest to carry out an extensive numerical study of the computational complexity of 
the schemes presented in this paper. This will be the subject of future research. 
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5 Appendix 

5.1 Proof of Lemma [3] 

We first suppose that p = 1. According to Theorem 5.2 page 72 of Milstein [1995], it suffices to check 
that there exists a positive constant C independent of N such that 
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Thanks to Ito's formula and to assumption (?4U: we have that 
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Using assumptions an d ("ft©, we also have Vp > 1 
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This implies both the second and the third inequality of (|33|) . This estimation is also sufficient to extend 
the result of Milstein [1995] to the L 2p norm and conclude the proof. 



5.2 Proof of Lemma [8] 

One can easily check that (Y t )o<t<T is a Gaussian process which has the same distribution law as the 
process (y e~ Kt + 9(1 - er Kt ) + ^£w#«t-i) Q < t <T. So, 



E 



J Po<t<T \Y t \ 1 + a 



< CE ^ e CSU P0<t<T \ W a^t _ 1 \ 1 + " 2 



Since sup 0<f<e 2 K T_ 1 \ Wt\ — (sup 0<t<e 3„T_ 1 Wt) V (— inf 0<t<e 2 K T_ 1 Wt), we deduce from the symmetry prop- 
erty of the Brownian motion that 

E (>isup < t < T |r t | 1 +^ < CE(e clsup °<t<° 2 « T -i Wtll+C2 + e c ' |inf °<*<= 2 ' sT -i H/ * |1+C2 ) 

< 2CE ^e c ' sup o<t<<! 2 ' tT ~i w t\ 1+ " 2 " 



The probability density function of sup 0<t<T Wt is equal to y i— > J 2t t{ y>Q } (see for example problem 
8.2 p. 96 of Karatzas and Shreve [1991]) which permits to conclude. 



5.3 Proof of Lemma 

The first point is an obvious consequence of the Feynman-Kac theorem. In order to prove the second 
one, let us first check the following result : 



For any multi- index (3 £ IN 3 such that f}\ < 6, 3 Cp, Kp > and pp £ N such that 
V(»,m,») e V T , \d plx (y,m,v)\ < Cpe~ K ^ (1 + \x\^) 



(34) 
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Indeed, using Leibniz's formula, one can show that dpj x (y, m, v) can be written as a weighted sum of terms 
of the form 



(x - log(ao) + pF(y ) - pF(y) - m) k * 
Ck = 7—1 exp 



2(1 -p> 



where fc = (ki,k2-,k^) belongs to a Unit set 7g C N 3 and (ai)o<i<fc 3 are constants taking value in {0,1}. 
Using assumption (T ST^ and ff/ TPfj) and Young's inequality, we show that 3Ck,Kk > and € N such 
that |0fe| < Cfee _ifkx2 (l + |a;| Pfc ) which yields the desired result. 

Now, let us fix a S N 3 , 1 € N such that 2/ + |a| < 6 and (t,y,m,v) G [0,T] x V t . Thanks to PDE (j24|) . 
d l t d a u x (t,y,m,v) = (—l) l d a C l u x (t,y,m,v). One can check that the right hand side is equal to a weighted 
sum of terms of the form dpiU x (t, y, m, v) x np a (b, a, f, h) where j3\ £ N 3 is multi-index belonging to a finite 
set J* ; , /?2 is a suffix belonging to a finite set I 2 l and (6, <J, /, /i) is a product of terms involving the 
functions 6, er, /, ft and their derivatives up to order 4. 

On the first hand, assumptions (T ST2^ and (7 fl3)) yield that 3cf > and qi jOC G N such that 



V/?2 € 72,1. kfc(M,/,&)l < cL(l + \y\«~). 



(35) 



On the other hand, by inverting expectation and differentiations, we see that d^Uxit, y, m, v) is equal 
to the expectation of a product between derivatives of the flow (y,m,v) — > (Yr-t,mT-t,'VT-tP v ' m ' v ' an d 
derivatives of the function ^ x evaluated at {YT~t,'mT-t,VT~tY y,m ' v ^ € T>t- Using result (|34l) and the fact 
that, under assumptions (7ifT2|) and (7i ll3|) . the derivatives of the flow satisfy a system of SDEs with Lipschitz 
continuous coefficients (see for example Kunita [1984]) we show that 3cj a , Ki^ a > and pi >a € N such that 



Vft e III, \d^u x (t,y,m,v)\ < cl a e- K ^(l + \x\^). 
Gathering ([35]) and (f36|) enables us to conclude. 



(36) 



5.4 Proof of Lemma _ 

Making the link between ODEs and SDEs (see Doss [1977]), one can check that (Y^, . . . , Y t ) has the 
same distribution law as (Y2ti > • ■ • !^2tjv) where (Yt)te[o,2T) 1S solution of the following inhomogeneous SDE 
% = yo + J*b(s,¥ s )ds + f* a(s,¥ s )dW s with, 
V(«,y)G [0,2T]x R, 



b{s,y) 



N-l r 



Ky)-\°°'{y) if se [j 



(y) 



fc=0 

otherwise 



(4fc+l)T (4/c + 3)T 



27V 



27V 



and 



N-l 



<?(s,y) = 



o ifse \J 

fc=0 

cr(y) otherwise 



(4fc + l)T (4fc + 3)T 



2iV 



2iV 



Since these coefficient have a uniform in time linear growth in the spatial variable, one easily concludes. 
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